Equation-free dynamic renormalization in a glassy compaction model 
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Combining dynamic renormalization with equation-free computational tools, we study the appar- 
ently asymptotically self-similar evolution of void distribution dynamics in the diffusion-deposition 
problem proposed by Stinchcombe and Depken [Phys. Rev. Lett. 88, 125701 (2002)]. We illustrate 
fixed point and dynamic approaches, forward as well as backward in time. 
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The study of glasses is an important topic under in- 
tense investigation for several decades now (see e.g., the 
review of Experimental studies in granular com- 

paction (such as the ones in have offered consid- 
erable insights on the behavior and temporal dynamics 
of glass forming systems. From the theoretical point of 
view, a variety of modeling/computational approaches 
have been used to understand better the key observables 
and their evolution for this slow relaxational type of dy- 
namics. While direct (molecular dynamics) simulations 
are widely used, perhaps more popular are kinetic ap- 
proaches including e.g., the Fredrickson- Anderson model 
[sj and more recent variants thereof 0. Another com- 
mon approach involves the use of mode coupling theory 
examining the time evolution of the decay of density 
fluctuations to separate liquid from glassy (non-ergodic 
due to structural arrest) dynamics. More recently, the 
energy landscape of glass forming systems and the role 
of its "topography" has become the focus of numerous 
works 6] . 

In a previous paper we proposed a simple com- 
paction "thought experiment" for hard spheres: the in- 
sertion of a hard sphere in a gas of hard spheres (accepted 
when the sphere does not overlap with previously exist- 
ing ones). Combining simple thermodynamic arguments 
with equilibrium distribution results we argued that the 
evolution of the hard sphere density should be logarith- 
mic in time (as the maximal density is approached). 

In a recent paper Stinchcombe and Depken pre- 
sented an interesting diffusion-deposition model exhibit- 
ing "glassy" compaction kinetics (see also :9j for a similar 
model example). The model provides a useful paradigm 
for testing the hypothesis of self-similar evolution of the 
system statistics -in particular, of the density and the 
void distribution function- based on direct simulations. 

The techniques that we will use to test this hy- 
pothesis combine template-based dynamic renormaliza- 
tion with the so-called "equation-free" approach to 
complex/multi-scale system modeling Dynamic 
renormalization has been developed in the context of 
partial differential equations with self-similar solutions 
(e.g. blowups in finite time) 0|. The equation-free 



approach targets systems for which a fine scale, atom- 
istic/stochastic simulator is available, but for which no 
closed form coarse-grained, macroscopic evolution equa- 
tion has been derived. The main idea is to substitute 
function and derivative evaluations of the unavailable 
coarse-grained evolution equation with short bursts of ap- 
propriately initialized fine scale computation. The quan- 
tities required for traditional continuum numerical anal- 
ysis are thus estimated on demand from brief computa- 
tional experimentation with the fine scale solver. We will 
demonstrate equation-free fixed point approaches to com- 
puting the self-similar shape and dynamics of the void 
cumulative distribution function. 

We start with a brief description of the model and 
key (for our purposes) results of Q. We then perform 
equation-free dynamic renormalization computations for 
this problem. The coarse-grained description of the sys- 
tem evolution backward in time is examined through re- 
verse coarse projective integration. We conclude with a 
brief discussion of potential extensions of the approach. 

The model considered by Stinchcombe and Depken 
consists of unit-size grains interacting through hard-core 
potentials and performing a Monte Carlo random walk. 
While some of their results apply in any dimension, in 
this paper we will work in one spatial dimension; this is 
for convenience-our approach is not limited to Id. While 
the grains diffuse on the line, as soon as a sufficiently 
large void is formed, it is instantaneously filled by an 
additional grain. We will work with a system of finite 
size and periodic boundary conditions; the system size L 
(here of the order of 10^) is large compared to the grain 
size, which is taken to be 1. If p is the system density, 
(5 = (1 — p)/ p, is the average void size. One of the main 
findings in 8] is that the void density e ^ \ — p goes to 
zero as an inverse logarithm 
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A similar asymptotic behavior 



asymptotically for large t. 
follows for 5. 

Density provides a cumulative measure of the system 
evolution (it is the zeroth moment of the particle distribu- 
tion function) . We set out to study in detail the evolution 
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FIG. 1: Dynamic evolution of the model for different initial- 
izations. Dotted (resp. dash-dotted) lines: starting with the 
self-similar (resp. uniform) void CDF. Solid (resp. dashed) 
lines: evolution replotted in terms of an appropriate "initial 
time" to; notice the collapse of the curves on a straight line, 
see text. The inset shows void CDFs at different time/density 
instances (color-coded as in the main Figure) ; they all collapse 
upon rescaling to a single CDF (smaller inset). 
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FIG. 2: Schematic of coarse dynamic renormalization. 



of the void cumulative distribution function, C'DF{x,t). 
This function is defined between and 1 in the void size 
X, since at no time do we get voids larger than 1 (the mo- 
ment they appear, they get filled). The dotted lines (and 
the dash-dot line) in Fig. ^show the evolution of 6{t) 
for various initial void distributions (the choice of initial 
void distributions for these transients will be discussed 
below). The stochastic simulation uses 100,000 particles; 
sequential update is used, with a maximum spatial step 
size of 0.5. The semilog plot clearly indicates the asymp- 
totic logarithmic regime in time. The inset shows the 
shape of the void CDF for various simulations and vari- 
ous instances in time, within the logarithmic regime; all 
shapes when appropriately rescaled with the average void 
size collapse on a single curve (smaller inset). These two 
observations (logarithmic time dependence, and collapse 
of the distributions upon rescaling) clearly suggest the 
possibility of self-similar evolution dynamics. 

In particular, the graph of the time evolution of the 
densities suggests that S becomes eventually proportional 
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Replotting the data in view of the above relation col- 
lapses the Fig. n curves onto a single line; the role of to 
will be discussed below. This implies that the dynamical 
evolution of 5 eventually follows the differential equation: 
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(and hence a similar equation can easily be obtained for 
e). The dynamics of ^ share the asymptotic behavior 
ofEq. 111). 

We now test the self-similarity hypothesis for the CDF 
evolution. If the dynamics of the CDF are self-similar 
and stable, all initial distributions will asymptotically 
approach the same self-similar shape (modulo rescaling) 
while they densify and slow down. Our goal is to find 
the self-similar shape at a convenient scale - a relatively 
low density, when the evolution dynamics are still fast. 
Our main assumption is that the void CDF is a good 
observable for the system dynamics - i.e., that a macro- 
scopic evolution equation (possibly averaged over sev- 
eral experiments) conceptually exists for its dynamics, 
even though it is not available in closed form. Starting 
with a given void CDFq with 6 — So we lift to an en- 
semble of microscopic realizations of it - that is, grain 
configurations possessing the given void CDF. This 
"lifting" is accomplished here by randomly selecting the 
voids from CDFq as grains are "put down" on the line; 
Sq indirectly selects the scale at which we will perform 
our computation. We then evolve each of these real- 
izations based on true system stochastic dynamics for 
a time horizon r. Finally, we obtain the ensemble av- 
eraged void CDF{x,t) and its S{t) = 5i. This is the 
restriction step of equation-free computation: evaluating 
the macroscopic observables of detailed, fine scale com- 
putations. We now rescale our macroscopic observable 
(the CDF{x,t)), using the ensemble average Si, to the 
original ,5o: CDF{Sox/Si,t) ^ CDF{x,t). (Clearly this 
requires the largest void in the rescaled CDF to be less 
than 1 - a condition that one may expect to prevail at 
high enough densities) . We then discard the simulations 
we performed, and start a new ensemble of simulations at 
the original density, but with the more "mature" CDF. 
The map from current void CDF to future void CDF 
is the "coarse timestepper" of the unavailable equation 
for the macroscopic observable evolution; the composi- 
tion of this map with the rescaling step constitutes the 
"renormalized coarse timestepper" , the main tool of our 
dynamic renormalization procedure, schematically sum- 
marized in Fig. 121 Given this map, several approaches 
to the computation of the long-term coarse self-similar 



3 



dynamics exist. The simplest is successive substitution 
- we repeat coarse renormalized time-stepping again and 
again, and observe the approach of the void CDF to its 
self-similar shape at the scale we chose (parametrized by 
60) ■ The fixed point of this procedure lies on the group 
orbit of scale invariance for the CDF dynamics; it was se- 
lected through a "pinning" or "template" condition (our 
choice of Sq). This dynamic renormalization iteration, 
with T = 4 is shown in Fig. starting from a uni- 
form void distribution with So = 0.2; the inset shows the 
third iteration before (dotted line) and after (solid line) 
rescaling. 
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FIG. 3: a. Successive evolution of an initially uniform CDF 
profile with renormalization: the inset shows the profile at 
the n=3 iteration before (dashed) and after (solid) rescaling. 
b. Newton-Krylov-GMRES fixed point renormalization cal- 
culation, starting from uniform void CDF; the inset shows the 
decrease of the norm of the residual vector at each Newton 
step. 

Successive substitution will converge to stable self- 
similar solutions. It is also, however, possible, to use 
fixed point algorithms to converge to the self-similar void 
CDF at So- In effect, we wish to solve a functional equa- 
tion for the fixed point of the renormalized timestepper: 

CDF{x,t) -CDF{x,0) ^0. (4) 

Fixed point algorithms (like Newton's method) require 
the repeated solution of systems of linear equations in- 
volving the Jacobian matrix of some discretization of 
Eq. Since we have no closed form equations for 

the void CDF evolution, this Jacobian is unavailable; 
yet equation-free methods of iterative linear algebra (like 
GMRES 0) allow the solution of the problem through 
a sequence of matrix-vector product estimations. In our 
case these estimations are obtained through the "lift-run- 
restrict-rescale" protocol performed at appropriately se- 
lected nearby initial distributions (for details, see [ll|). 
Short bursts of stochastic simulation from nearby initial 
void distributions allows us to estimate the action of the 
Jacobian on selected perturbations, and hence the solu- 
tion of linear equations and ultimately, of the nonlinear 



fixed point problem. Fig. |3Jd shows the iterates of such a 
matrix-free fixed point computation for Eq. Q through 
a Newton-Krylov GMRES algorithm applied to a 100- 
point uniform finite difference discretization of the void 
CDF in x; our initial guess was a uniform void CDF 
{So = 0.2), and we converge to the self-similar shape 
within a few iterations; the inset shows the evolution 
of the norm of the residual vector with iteration num- 
ber. Lifting from this coarse description (100 numbers) 
to realizations of the CDF (200,000 particles) involved 
linear interpolation of the CDF. Relatively short runs 
(r = 1) were used to construct the coarse renormalized 
timestepper in Fig. |3t); its fixed point is independent of 
the reporting horizon r. 

Wrapping traditional numerical algorithms around 
coarse timesteppers enables several computational tasks, 
of which the fixed point computation above is only one 
example. Another interesting example is the so-called 
coarse projective integration; here, short runs of the 
stochastic simulation are used to estimate the local time 
derivatives of the macroscopic observables (e.g. of the 
void CDF); these estimates are used to "project" the 
void CDF "far" into the future. The simplest such algo- 
rithm is projective forward Euler; see [llj for multi-step 
and implicit coarse projective integration. One does not 
evolve the microscopic problem; one evolves the coarse- 
grained closure of the microscopic problem, which is pos- 
tulated to exist, but is not explicitly available. Here we 
apply coarse projective integration backward in time (as 
discussed in jlj this is appropriate for systems with large 
time scale separation between fast stable dynamics and 
slow dynamics). Starting with the self-similar CDF at 
some (relatively large) density (see Fig. 0]and its insets) 
we lift to microscopic realizations, evolve briefly forward 
in time, and estimate the rate of void CDF evolution. 

This rate (estimated from the short pink run) is used 
to project the void CDF for a (longer) interval backward 
in time (in blue). The forward simulations are then dis- 
carded, and new simulations are initialized at the "ear- 
lier" void CDF. Repeating the procedure clearly shows 
an acceleration in backward time as the density becomes 
lower; the probability of large voids becomes increasingly 
larger, the backward dynamics evolves faster and faster 
and their fluctuations intensify, making rate estimation 
from the short forward runs difficult. Fig. ^seems to sug- 
gest a "reverse finite time" rate singularity. This is only a 
visual suggestion, however; as the density becomes lower, 
we do not expect an evolution equation for the void CDF 
to provide a good model. More detailed physical mod- 
eling of the void-filling process will be required for the 
study of this backward dynamics at lower densities to 
become meaningful. 

It is this "visual suggestion" of a backward explosion 
of the density evolution rate that suggested fitting the 
data as a function of t + to', we find that, if we initialize 
with the self-similar CDF (at whatever density), evolu- 
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FIG. 4: Reverse coarse projective integration: the dotted pink 
line is the forward estimation of the rate; the solid blue line is 
the projection backward in time. The upper inset shows (in 
the same color code as the large dots in the figure) the CDF 
profile at the beginning and the end of the overall run. 

tion data can be fitted almost perfectly with a straight 
line for the appropriate tg in Fig. ^ Based on Eq. 
the value of the term B from all of our trajectories sug- 
gests that at time t + to = 1 the self-similar CDF density 
is ^ 0.60; that is, tg — 1 is the time that it takes for a 
simulation initialized with the self-similar void CDF at 
density 0.60 to evolve to the initial density of our com- 
putational experiments (our i = 0). When the initial 
condition is different than the self-similar CDF (e.g. if 
it is a uniform distribution, see the dash-dot line in Fig. 
^ finding a good t^ does not collapse the entire transient 
on a straight line - we see, however, that the solution 
asymptotically approaches the self-similar one (dashed 
line). 

In summary, in this paper we implemented equation- 
free, coarse-grained dynamic rcnormalization (simulation 
forward and backward in time, as well as fixed point 
computations) to study the evolution of the void CDF 
for a model glassy compaction problem; we found this 
evolution to be governed by apparent asymptotic self- 
similarity over the range of densities we could reliably 
compute. The procedure is, in principle, capable of con- 
verging to both stable and unstable self similar solutions, 
find similarity exponents when they exist, as well as to 
quantify the fixed point stability (by using matrix-free it- 
erative eigensolvers like Arnoldi) . In a continuation / bi- 
furcation context the algorithms can be used to track self- 
similar solutions in parameter space, detect their losses 
of stability and bifurcations; of particular interest is the 
onsets in parameter space, of self-similar dynamics, which 
will appear in our formulation as a fixed point bifurcation 
Finally, coarse projective integration (forward or 
reverse) can be performed in physical space (void CDF 



evolution in time) or in renormalized space (renormal- 
ized void CDF evolution in logarithmic time). We be- 
lieve that the bridging of continuum numerical techniques 
with microscopic simulations we illustrated here may be 
helpful in the coarse grained study of glassy dynamics 
through atomistic/stochastic models, especially when the 
dynamics depend on parameters of the microscopic rules. 
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